High resolution radon transform seismic traces processing

ABSTRACT

The present invention comprises a non-iterative method of processing seismic traces. A constrained High Resolution Radon decomposition is performed at various frequencies in which the Radon decomposition at a given frequency is constrained as a function of the Radon decomposition at at least a lower frequency. It is emphasized that this abstract is provided to comply with the rules requiring an abstract which will allow a searcher or other reader to quickly ascertain the subject matter of the technical disclosure. It is submitted with the understanding that it will not be used to interpret or limit the scope of meaning of the claims.

This application claims the benefit of Provisional application Ser. No60/132,211, filed May 3, 1999.

BACKGROUND OF THE INVENTION

1. Field of the invention

This invention relates to an improved high resolution Radon transformfor use in analysing geophysical data.

2. Description of the Prior Art

The widespread use of multiple removal in the parabolic Radon domain(Hampson, 1986, Kostov, 1990) is related to its effectiveness andefficiency for most situations. However, when applied under severeconditions reduced spatial aperture (FIGS. 1a-f) and/or coarse offsetsampling (FIGS. 2a-c), one may observe a poor focusing of the events inthe parabolic Radon domain, combined with severe allasing artifacts. Asa consequence, the multiple model loses some multiple energy and includesignificant primary energy. Once the multiple model is substracted fromthe input data, this leads to poor multiple removal and deteriorateprimaries (see FIGS. 3a-d for illustration).

De-aliased Hiqh-Resolution (DHR) Radon Transform

Finite spatial aperture limits the resolution of the Radon transform,while finite spatial sampling introduces aliasing artifacts. To overcomethese limitations one has to constrain the parabolic decomposition ofthe data. This issue was first investigated by Thorson and Claerbout,1985. More recently Sacchi and Ulrych, 1995, Hugonner and Canadas, 1997,and Cary, 1998 have developed high-resolution Radon transforms (in thefrequency-space or time-space domain). These constrain the Radon spectrato be sparse in q and t, using a re-welghted iterative approach.

However, such an iterative approach presents the major drawback of beingvery much time and computer consuming.

SUMMARY OF THE INVENTION

The present invention is directed toward a novel method that has theadvantage of not being an iterative process.

Going back to the previous synthetic example (FIGS. 3a-c), the followingobservations can be made:

1. Due to the curvature range involved in this example and the parabolicsampling rate chosen, 550 parabolas are used to perform the parabolicdecomposition of 15 traces. That is an under-determined least-squaresproblem.

2. Among the 550 parabolas, only 4 are actually needed to properlydecompose the data. Constraining the parabolic decomposition of the dataonto these four parabolas will lead to well-focused parabolic Radonspectra. One then has to solve a constrained underdeterminedleast-squares problem.

3. At low frequencies, the steering vectors used for the parabolicdecomposition do not suffer from aliasing (FIG. 2b). As a consequencethe parabolic decomposition at low frequencies can be used to guide theparabolic decomposition at higher frequencies.

To handle this constrained under-determined least-squares problem, theinvention proposes a data driven constrained Radon decomposition. TheRadon decomposition at a given frequency ω_(k) is constrained around themost significant spectral components highlighted at the previousfrequency ω_(k−1). This non-iterative, gradual way (from low frequenciesto high frequencies) to build the constrain enables to enhance theresolution of the Radon spectra. This algorithm enables to go beyond thecommonly admitted sampling and aperture limitations. The task of theproposed method is simplified when the data to decompose are solelycomposed of a small amount of Radon components. On actual data thisapproach is therefor more effective using sliding temporal and spatialwindows.

More generally, the invention proposes a method of performing aprocessing on seismic traces comprising the step of performing aconstrained High Resolution Radon decomposition at various frequencies,wherein the Radon decomposition at a given frequency is constrained as afunction of the Radon decomposition at at least a lower frequency.

More particularly, it proposes a method wherein a Radon decomposition issuccessively performed at various sparse frequencies, from the lowerfrequency to the higher, the Radon decomposition at a given frequencybeing constrained as a function of the Radon decomposition at theprevious frequency.

Other features and advantages of the invention will be furtherunderstood in view of the following description.

DESCRIPTION OF THE DRAWINGS

The patent or application file contains at least one drawing executed incolor. Copies of this patent or patent application publication withcolor drawing(s) will be provided by the Office upon request and paymentof the necessary fee.

FIGS. 1a-b depicts the effects of the aperture limitation on thediscrimination of three parabolic events.

FIGS. 1e-f depicts the horizontal event resolution in the parabolicRadon domain (1 d, 1 f) as a function of the spatial aperture (1 c, 1e). The origin of the ‘bufferfly’ looking artifact in (1 f) can be foundin (Kabir and Marfurt, 1999).

FIG. (2 a) depicts a parabolic event for 0<f<90 Hz—same parabolic eventat 5 Hz (2 b) and at 70 Hz (2 c). Due to the finite spatial sampling,spatial aliasing patterns appear, leading to a non-unique parabolicdecomposition of this monochromatic wave-field. The parabolicdecomposition must be constrained by the parabolic decomposition atlower frequencies.

FIG. (3 a) depicts a simulated NMO-corrected CMP gather. (3 b) Multiplemodel obtained using regular parabolic Radon transform: we observe thepresence of primary energy. (3 c) Estimated primaries ((3 a)-(3 b)):note the deterioration of the primary energy and the remaining multipleenergy. (3 d) Regular parabolic Radon spectra.

FIG. (4 a) depicts the same simulated NMO-corrected CMP gather. (4 b)Multiple model obtained using DHR parabolic Radon transform. (4 c)Estimated primaries ((4 a)-(4 b)), note the perfect multiple removalwithout damaging the primaries. (4 d) DHR parabolic spectra.

FIG. (5 a) depicts an NMO-corrected CMP gather from deep offshoreexploration. (5 b) Multiple model estimated from a windowed portion ofthe CMP gather. (5 c) estimated primaries ((5 a)-(5 b)). (5 d) Estimatedprimaries after moving the sliding window all over the CMP gather.

FIG. 6 is a flowchart of an exemplary method of the present invention.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

It is now described how De-aliased, High-Resolution (DHR) Radon spectramay be constructed with a direct approach avoiding the disadvantagesassociated with the iterative high-resolution Radon transform of theprior art.

The Radon transform (linear of parabolic) allows to decompose a complex,signal as a sum of elementary signals with simpler spatial behaviors(linear or parabolic). Once the signal is decomposed in the Radon space,some of its spectral components can be filtered out in order to removesome unwanted coherent (linear or parabolic) components of the signal.The effectiveness of the signal filtering in the Radon domain relies onthe quality of the Radon transform: its ability to uniquely decomposethe signal (Radon spectra). The uniqueness of the Radon decomposition isaltered when the number of spectral components involved in the Radontransform exceeds the number of signal components. As a consequence ofthis non-uniqueness aliasing patterns appear in the Radon spectra. Thisaliasing phenomena is illustrated with FIG. (3 d), which represents theresult of the decomposition on 15 traces NMO-corrected (FIG. 3a) over550 parabolas. The nine parabolas contained in the input signal have notbeen uniquely decomposed, leading to unfocused Radon spectra.

In order to enhance the resolution of the Radon spectra one has toconstrain the parabolic decomposition of the signal onto its mostsignificant spectral. components. How to define and apply theseconstrains is now described.

For efficiency reasons the Radon decomposition is often performed in thefrequency domain. Given a frequency component ω_(k), the followingequality relates the recorded data to its spectral components:

{overscore (d)}(ω_(k))=G(ω_(k)){overscore (q)}(ω_(k))  (1a)

with ω_(k), the k^(th) frequency component, {overscore (d)}(ω_(k)) thedata vector to decompose (M: number of available data samples),

{overscore (d)}(ω_(k))=(d ₁(ω_(k)), . . . , d _(M)(ω_(k)))^(T)  (1b)

{overscore (q)}(ω_(k)), the spectral vector containing the N Radonspectral components over which the data vector is decomposed,

{overscore (q)}(ω_(k))=(q ₁(ω_(k)), . . . , q _(N)(ω_(k)))^(T)  (1c)

G(ω_(k)) the (M, N) complex matrix made of the Radon steering vectors,

G _(m,n)(ω_(k))=e ^(Jω) _(k) ^(q) _(n) ^(x) _(m) (Linear Radontransform)  (1d)

G _(m,n)(ω_(k))=e ^(jω) _(k) ^(q) _(n) ^(x) _(m) ^(₁) (Parabolic Radontransform)  (1e)

(x_(m): spatial co-ordinate at which d_(m)(ω_(k)) is recorded).

The Radon decomposition of the data is usually obtained using theover-determined Least-squares solution to (1a),

{overscore (q)}(ω_(k))=G ^(H)(ω_(k))G(ω_(k))+εI)⁻¹ G^(H)(ω_(k)){overscore (d)}(ω_(k))  (2a)

with ε a pre-whitening factor to avoid numerical instabilities and I the(N, N) identity matrix.

In order to constrain the Radon decomposition onto the most significantspectral components of the data, it is proposed to compute the Radondecomposition using the constrained under-determined Least-squaressolution to (1a),

{overscore (q)}=W(ω_(k))G ^(H)(G ^(H) W(ω_(k))G+εI)⁻¹ G ^(H) {overscore(d)}  (2b)

with W(ω_(k)) (N, N) a real diagonal positive definite constrain matrix,that will focus the Radon decomposition around the most significantspectral components of the data.

How to fill the W diagonal matrix will now be described.

Assuming, non-disperse signals, the Radon amplitude spectra's has somesort of continuity from one frequency to another one. This observationis used to constrain the Radon spectra at frequency ω_(k) with the Radonspectra derived at the previous frequency ω_(k−1). Therefore theconstrain matrix W(ω_(k)) reads,

W _(i,l)(ω_(k))=∥q _(i)(ω_(k−1))∥i=l, . . . , N  (3)

It will be readily understood that with such a constrain matrixW(ω_(k)), equation (2b) can easily be solved using for example thealgorithm proposed by Sacchi and Ulrych.

After the processing of the Radon decomposition, the data are filteredin the Radon space to substract the multiples and an inverse parabolicor linear Radon transform is performed to obtain the estimatedprimaries.

This non-iterative, gradual way (from low frequencies to highfrequencies) to update the constrain matrix enables to enhance theresolution of the Radon decomposition and to avoid a large amount ofaliasing artifacts.

In particular, as clearly shown by the example displayed on FIGS. 4a-4 dwithout a prioi information on the curvature of the multiples, thisnon-iterative process focuses the parabolic decomposition of the dataonto its most significant spectral components. The application of thismethod to the previous synthetic data example leads to remarkableresults (FIGS. 4a-d), including a sparse parabolic decomposition of thedata along the q and τ axis, leading to perfect multiple removal with nodamage to the primaries.

Application to Real Data

The task of the proposed method is simplified when the data are solelycomposed of a small number of parabolas. On actual data this approach ismore effective using sliding temporal (200 ms) and spatial windows (18traces) as displayed in FIGS. 5a-d. This example nicely illustrates theability of the algorithm to separate, on a limited number of traces,primaries from multiples with large move-out The algorithm has gonebeyond the usual sampling and aperture limitations.

The proposed non-iterative De-aliased, High Resolution Radon transformof the present invention provides an alternative to the traditionalRadon transform when one has to handle severe circumstances; smallspatial aperture, insufficient spatial sampling, large or small move-outdifference between primaries and multiples. Working on limited spatialand temporal windows, the wave-field is readily decomposed into its mainparabolic components using the present invention.

REFERENCES

Cary, P., 1998. The simplest discrete Radon transform. ExtendedAbstracts, Vol. II, p. 1999-2002.

Hampson, D., 1986. Inverse, velocity stacking for multiple elimination.J. Can. SEG, 22, p. 44-45.

Hugonnet, P. and Canadas, G., 1997. Regridding of irregular data using3D Radon Decompositions: SEG Extended Abstracts, Vol. II, p. 1111-1114.

Kabir, M. M. N. and Marfurt, K. J., 1999. Toward true amplitude multipleremoval. The leading Edge, Vol. 18n N1, p. 66-73.

Kostov, C., 1990. Toeplitz structure in Slant-Stack Inversion: SEGAbstract Vol. II, p. 1647-1650.

Sacchi, M. D. and Ulrych, T. J., 1995. High resolution velocity gathersand offset-space reconstruction: Geophysics, 60, 1169-1177.

Sacchi, M. D., SEG Expanded Abstracts 1999, Fast High resolutionparabolic Radon Transform.

Spitz, S., 1991. Seismic trace interpolation in the F-X domain.Geophysics, Vol. 56, N6, p. 785-794.

Tarabtola, A., 1987. Inverse Problem Theory Elsevier.

Thorson, J. R. and Claerbout, J. F., 1985. Velocity-stack andslant-stack stochastic inversion: Geophysics, 50, p. 2727-2741.

The foregoing disclosure and description of the invention areillustrative and explanatory. Various changes in the details of theillustrative embodiments may be made without departing from the spiritof the invention.

What is claimed is:
 1. A method of performing a processing on seismictraces comprising performing a constrained High Resolution Radondecomposition at various frequencies, wherein the Radon decomposition ata given frequency is constrained as a function of the Radondecomposition at at least a lower frequency.
 2. A method according toclaim 1, wherein a Radon decomposition is successively performed atvarious sparse frequencies, from the lower frequency to the higher, theRadon decomposition at a given frequency being constrained as a functionof the Radon decomposition at the previous frequency.
 3. A methodaccording to claim 1, wherein the constrained Radon decomposition isprocessed using a constrained over-determined Least-squares solution tothe equation: {overscore (d)}(ω_(k))=G(ω_(k)){overscore (q)}(ω_(k)) withω_(i), the i^(th) frequency component, {overscore (d)}(ω_(i)) the datavector to decompose, {overscore (q)}(ω_(k)), the spectral vectorcontaining the Radon spectral components over which the data vector isdecomposed.
 4. A method according to claim 2, wherein the constrainedRadon decomposition is processed through the determining of theLeast-squares solution to the equation: {overscore (q)}=W(ω_(k))G ^(H)(G^(H) W(ω_(k))G+εI)⁻¹ G ^(H) {overscore (d)} with W(ω_(k)) a realdiagonal positive definite constrain matrix which components arefunction of at least a frequency lower than ω_(i).
 5. A method accordingto claim 4, wherein the constrain matrix W(ω_(k)) reads, W_(i,i)(ω_(k))=∥q _(i)(ω_(k−1))∥.
 6. A method according to any of thepreceding claims, wherein the Radon decomposition is filtered in theRadon domain.
 7. A method of performing a processing on seismic traceswherein the method according to any of claims 1 through 5 is applied onsliding temporal and spatial windows.